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ABSr  ACT 


A  three  dimensional  hyperbolic  differential  equation 
based  on  finite  correlated  particle  velocities  Is  derived 
which  is  appropriate  to  modeling  anisotropic  turbulent 
diffusion  in  the  atmosphere.  Cauchy  initial  data,  the 
mean  wind,  the  Reynolds  stress  tensor,  and  a  typical 
frequency  of  pulsation  are  required  fcr  complete  solution. 
The  outlines  of  plumes  and  puffs  may  be  obtained  with  only 
knowledge  of  the  Reynolds  stress  tensor  and  mean  wind 
velocity.  The  classical  parabolic  diffusion  equations  are 
a  limiting  form  of  this  hyperbolic  model. 
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I.  INTRODUCTION  TO  ATMOSPHERIC  TURBULENT  DIFFUSION 


This  phase  space  derivation  of  a  hyperbolic  turbulent 
diffusion  equation  Is  motivated  by  several  less  general 
random  walk  and  hyperbolic  diffusion  models  [6,  8,  13, 

15,  16,  18,  19*  263  which  appear  to  more  reasonably  re¬ 
present  the  physical  processes  involved  In  atmospheric 
diffusion  than  any  of  the  classical  models.  The  basis 
for  the  classical  models  of  turbulent  diffusion,  almost 
without  exception,  is  one  form  or  another  of  the  parabolic 
diffusion  equation.  Even  the  seml-empirical  Gaussian 
plume  equation,  the  foundation  for  much  recent  diffusion 
studies,  may  be  derived  from  a  parabolic  diffusion 
equation  [3],  Recent  reviews  of  the  state  of  diffusion 
modeling  In  the  lower  atmosphere  express  a  dissatis¬ 
faction  with  the  degree  of  generality  to  which  existing 
models  may  be  extended,  and  agree  that  no  comprehensive 
mechanism  of  turbulent  diffusion  has  been  modeled  [3, 

5,  1^].  It  Is  significant  that  most  recent  reviews 
concerned  with  diffusion  application  omit  serious 
consideration  of  hyperbolic  systems  such  as  those  of 
Monin  [19]  and  Goldstein  [13]  which  are  expressions  of 
particle  motions  correlated  as  in  the  case  of  atmospheric 
turbulence.  While  these  first  efforts  are  not  very 
general,  they  do  provide  certain  features  which  allow  a 
more  realistic  modeling  compatible  with  practical 
applications.  In  the  present  hyperbolic  model,  the 
primary  assumption,  that  diffusion  can  be  represented  by 
a  Markov  process  in  phase  space,  is  suggested  as  being 
a  less  restrictive  assumption  than  the  corresponding 
assumption  used  in  the  derivation  of  the  classical  para¬ 
bolic  diffusion  equation.  The  classical  parabolic 
diffusion  equation  can  be  derived  by  assuming  a  Markov 
process  in  configuration  space  C1*  10 ,  173 •  Further¬ 

more  the  classical  assumption  that  the  velocity  necessary 
for  the  formulation  of  a  continuity  equation  can  be 
derived  from  the  particle  density  gradient  presents  the 
paradoxical  result  that  densities  exist  everywhere  even 
a  short  time  after  a  release.  Thus  particle  paths  exist 
for  which  there  are  no  corresponding  velocities.  In  the 


V,  ♦**•*«►  t,  *»«*«*«*£ 


phase  space  model  the  particles  have  velocities  -  finite 
velocities.  The  possibility  exists  then  to  relate  the 
velocities  of  the  diffusing  particles  to  the  velocities 
of  the  turbulent  atmospheric  fluid. 
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II.  THEORETICAL  DERIVATION  OP  TURBULENT  DIFFUSION  EQUATION 


The  derivation  Is  based  on  two  premises.  The  first 
Is  that  turbulent  diffusion  can  be  described  by 
statistical  methods  relating  particle  statistics  to  the 
statistics  of  the  atmospheric  turbulence  field.  The 
second  is  that  the  mean  wind  profile  can  be  introduced 
into  the  turbulence  field. 

The  turbulence  field  representation  will  be  based 
on  the  "phase  space"  description  of  the  probability 
density  for  a  particle.  We  assume  that  the  diffusion 
can  be  represented  by  a  Markov  process  in  phase  space. 

The  phase  space  representation  of  a  particle  has  been 
investigated  by  Chandrasekhar  [l]f  Kramers  [17],  Obukhov 
[22],  Tchen  [29],  and  Davies  [8].  We  shall  start  similarly 
to  Chandrasekhar.  In  the  Markov  process  In  phase  space, 
the  probability  density  at  t+At,  W(t+At),  Is  derived 
from  the  probability  density  at  t,  W(t). 

+  -4  +■  4  +  -4  -f  4  4  4 

W(r,ujt+At)  -  //W(r-Arfu-Au,t)*(r-Ar,u-Au,Ar* Au)d( Au)d(  Ar) 


A  Is  the  transition  probability  density  in  phase  space,  r 
Is  tlje  vector  representation  of  the  configuration  space, 
and  u  is  the  vector  representation  of  the  velocity  space. 

We  now  relate  the  increments  of  time  and  configuration 
space  by 

♦  4- 

Ar  •  uAt •  (2) 


The  transition  probability  density  in  phase  space  can  now 
be  represented  by  the  transition  probability  density  in 
velocity  space  and  the  Dirac  delta  function. 

4  4  4  4  4  -4  -4 

*  ■  #(r-Ar,  u— Au, Au) f { Ar—uAt )  (3) 

The  Integration  over  the  configuration  space  Increment  is 
now  easily  performed. 

4  4  4  +«4  -4  -4  -4-4  -4  -4  -4 

W(r+uAt,u,t+At)  »/W(r,u-Au,t)  *{r,u-Au,Au)d( Au)  (4) 
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Both  sides  are  expanded  In  a  Taylor  series  about  the  point 
(r,  u,  t). 

W  ♦  «tAt  +  (u*  flOAt  ♦  0( At2 )  - 


♦  ...  } 


sw  ia2«  .  a2u 

/{W-E— —  AU. AU.  )2+  E  j— y-*AU  AU 
JUJ  1  2*Uj»  1  •u1*u^  1  J 

1<j  (5) 

{*-r|^-AU.4-i  i|JL(au  )2  +  E  -AU  AU  +  ...  }d( AU) 

i*“i  1  2  ^U^  1  i<j9ulauj  1  i 


For  the  moments  of  the  velocity  Increments  we  write: 

+«.  ^  ^ 

<AU,  >  ■  I  Au,  td(Au)  ;  <Au.  AU.  >  »  /4l.  Au  ed(AU) 

1  —  1  1  J  — 1  J 


(6) 


The  expansion  can  then  be  rewritten: 

,a2w 


V*  *  “•™lt  '  -  =5tSi<4V+  5  jliip<4ui2>-1Iwi3^Jli 

♦  Hfe  *  ,5.  531  i^>+5  <« 


♦  0(  ^U^AUj  AUJc>) 


Keeping  terms  to  first  order  In  At,  dividing  by  At,  and 
going  to  the  limit  of  At  *0  yields 


2 

W  eiMjW+vsMSw)  +  lim  o  ....  -  0 

t  At+O 


where 


t  m  11,8  <AU> 
a  *  A  t-^o  TT 


(8) 


(9) 
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Here  we  consider  a  as  a  mean  acceleration  with  position 

r,  and  velocity  3.  Thus  the  equation  of  continuity  In 
phase  space  Is 


Wt+u*?£W+$a*(5w)  -  o  (10) 

To  obtain  the  continuity  equation  in  configuration  space 
we  Integrate  over  velocity  space , 


Pt  +  v-(pU)  -  0  , 


(11) 


where 

and 


p  ■  ^Wdu  # 
p0  -  JJiWdu  . 


To  obtain  a  momentum  equation  we  multiply  the  continuity 
equation  in  phase  space  by  u  and  integrate  over  velocity 
space. 


(pU)t  +  $*(pttff)  -  pt  -  0  , 


(12) 


+» 

plJtJ  ■  /33Vdtf  , 

rnrnm 

pt  ■  /j|Wdu  . 

IJere  X  is  the  mean  exterior  force  per  unit  p  at  r,  t,  and 
U  and  $U  are  the  first  and  second^veloclty  moments 
respectively  of  the  particles  at  r,  t.  The  continuity 
equation  and  the  momentum  equation  now  form  a  hyperbolic 
system  of  first  order  partical  differential  equations  with 
the  dependent  variable  being  the  particle  density  in 
configuration  space.  This  system  of  partial  differential 
equations  is  the  diffusion  equation  in  configuration  space. 


where 

and 
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III.  INTRODUCTION  OP  MEAN  WIND  PROPILE  INTO  TURBULENCE 
FIELD 


To  Introduce  the  mean  wind  profile  into  the  turbulence 
field  we  represent  the  mean  velocity  of  the  diffusing 
particles  as  the  sum  of  the  mean  wind  and  a  macroscopic 
velocity  representative  of  the  turbulent  flux  with  reject 
to  the  mean  wind  field. 

0  -  V(x,  y,  z)  ♦  v(x,  y,  z,  t)  (13) 

Furthermore  we  assume  that  the  velocity  in  addition 
possesses  a  stochastic  term  which  may  be  represented 
a  trivariate  normal  distribution  in  the  components  of  u. 
This  distribution  is  constructed  such  as  to  require  that 
the  mean  outline  of  the  diffusing  particles  propagate 
with  a  velocity  dependent  on  the  stationary  wind  field 
independently  of  time.  A  distribution  satisfying  this 
requirement  is 

W  «{(2*r|M|  exp[-l/2X'M“1X]>*e(x,yfz,t).  (14) 


where  M  « 


W! 

V. 


11 

**12 

M13 

21 

U22 

I23 

•  “iT  K 

|i  -XX, 
12  12 

31 

W32 

u33 

etc. 


] M |  is  the  determinant  of  M  and  X'  is  the  row  matrix, 

X*  -  [xj,  *2»  X3I#  X  is  the  column  matrix  obtained  by 
transposing  X*.  We  have  xi«  u j—  Vj-  vlf  x2m  u2-  V2-  v2 , 

x3  ■  u3-  V3-  v3,  where_the  subscripts  1,  2,  3,  indicate 
the^coqponents  In  the  if  j,  k,  directions  respectively 
of  u,  V,  v.  It  is  clear  that  the  expected  value  03  is 


uu  - 


“21 

»*11 


W23I*J 

x£  v3  2k^ 

V33gk 

VV  +  Vv  + 


++  ■*-*  ■+-> 
vV  +  vv  +  00 


We  have  defined  33  as  the  tensor  of  variance s^and^co- 
variances  of  velocities  relative  to  'the  mean  V  +  v. 
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In  this  esse  an  assumption  that  (w  -foe)  i<vv»  is  in¬ 
dependent  of  time  gives  us  the  required  properties  of 
boundary  propagation.  We  define  pv  as  the  turbulent 
flux,  j,  and  p<^>  as  the  turbulent  energy  tensor.  The 
substitution  of  the  expressions  from  equations  (13)  and 

(15)  into  equations  (11)  and  &2)  yields 

pt  ♦  (pV+3)  -  0  (16) 

and 

(p^+S).  +v*[pVV^+v5+p<vv>]  -pX  -  0.  (17) 

Tf 

The  above  equations  constitute  a  hyperbolic  system 
of  first  order  partial  differential  equations.  These 
equations  describe  the  transfer  of  a  scalar  quantity, 
through  a  turbulent  fluid.  Particulate  diffusion  and 
heat  transfer  are  processes  for  which  the  equations  should 
hold. 

We  now  seek  to  relate^  and  «vv>  to  measurable 
quantities.  Once  X  and  <w>  are  given,  it  is  clear  that 
with  the  four  equations  (one  from  the  continuity  equation) 

(16) ,  and  the  three  from  the  vector  momentum  equation  (17), 
and  the  four  unknowns  (p  and  the  three  components  of  S), 
Cauchy  initial  data  is  sufficient  to  determine  the  hyper¬ 
bolic  system  of  equations  (16)  and  (17)  [2]. 

For  X  we  choose  an  acceleration  caused  by  a  resistive 
force  which  is  proportional  to  the  turbulent  flux  in  each 
direction  and  an  external  force  term  which  reflects  the 
contribution  of  gravity-buoyancy  effects  in  the  case  of 
heavy  particles. 

Ap  •  -  [  io(i*S)  +  J«(J*S)  +  kr(Jc*3)  3  ♦  «p  U8) 

The  resistive  force  coefficient  is  interpreted  as  the 
probability  per  unit  time  of  reversal  of  the  diffusing 
particle.  Nonin  [18],  Davies  [8],  and  Goldstein  [133 
Independently  derive  similar  expressions.  We  are  assum¬ 
ing  that  there  is  no  resistance  to  the  mean  flow,  v. 


Kasan^il^gnd  Mopin  report  Satisfactory  experimental 
verification ’of  a"si«4iiarf  ty  theory  application  using 
the  same  rfei stance  term  flp]  P2i)»  Monin  interprets 
the  resistance  term  as %  being:  equivalent  to  h(  typical 
frequency  of  pulsation  pjf\;the  turbulence,  Mbnin's* *  - 
approach  seems  consistent  hrttft  the  present  tievelbpment. 

The  values  of  this  term  mpy  be?  deduced  experimentally 
from  steady  state  evaluations -of* the  concentrations. 

For  the  turbulent  energy  tensor,  <$$>,  we  propose: 
that  the  Reynolds  stress  tenspxv,  he  used,  -We  propose 
that  the  second  velocity  moment  of1  the  diffusing  particles 
with  respect  to  the  mean  wind  field,  is  locally  pro¬ 
portional  to  that  of  the  turbulent  fluid.  At  present, 
for  simplicity  and  heuristic  purposes,  the  proportionality 
constant  will  be  assumed  as  unity.  Experiments  such  as 
smoke  puff  and  plume  outline  photography  may  determine 
more  appropriate  proportionality  coefficients.  It  must 
be  kept  in  mind  that  careful  consideration  of  the  type 
of  averaging  used  in  representing  <vv>  is  germane  to 
any  application.  The  <vv>  averaglng^must  he  compatible 
with  the  averaging  used  to  determine  V, 

Equation  (17)  then,  becomes 

+  «*[p$$+5$*$$+p<vv>3  + 


(19) 

Ia(i*8)  ♦  i«(j’S)  ♦  k  (k*3)  -  sp  "  Q. 

Formally  to  obtain  X  we  could  have  stalled  similarly 
to  the  Kramers  -  Chandrasekhar  approach  to  Brownian 
motion  in  a  field  of  force  [1],  The  increment  of  velocity 
which  the  particle  experiences  in  the  time  At  is 
expressed  &s  the  sum  of  a  term  due  to  the  external  field 
of  force,  Kat,  and  a  fluctuating  quantity  with  a  given 
la**  of  distribution,  S(4t) , 

a$  «  &  b  *  SUtK  UQ) 


The  classical  distribution  of  B(At)  Is  given  by 

(- j&(4t) | 2/4qat) 

^  1  « 
w(B[At])  ■  - - -  - 

( 4*qAt )  V*  (21) 


where 

q  «  BKT/m. 

Physically,  S(At)  represents  the  net  acceleration 
which  a  particle  experiences  In  time  at  under  the  In¬ 
fluence  of  molecular  scale  fluctuations.  We  shall 
neglect  It  compared  to  the  much  larger  scale  of  tur¬ 
bulence  fluctuation^.  Thus  $  also  reduce^  to  the  Dirac 
delta  function,  s(Au-KAt).  For  the  term  K  we  propose  a 
frictional  force  acting  on  the  diffusing  particle  pro¬ 
portional  to  the  turbulent  flux  velocity  and  external 
force  per  unit  p  due  to  buoyancy-gravity  effects. 

It  ■  -  [io(i'v)  ♦  J6(J*v)  ♦  ky(k*v)]  +  e 

(22) 

Substituting  these  values  into  equation  (4)  we  may 
formally  obtain  equations  (16)  and  (19)> 


IV.  NORMALIZATION 


In  any  application  of  the  hyperbolic  diffusion 
equation  it  must  be  remembered  that  p  is  a  probability 
density.  Hence  the  concentration  must  be  normalised. 

If  the  equation  were  parabolic,  a  single  tern  would  be 
sufficient  to  normalize  the  concentration.  However, 
for  the  hyperbolic  equation  it  is  necessary  to  perform 
two  integrations  to  normalize,  one  over  the  characteristic 
surfaces,  and  one  over  the  volume  contained  therein. 

Davies  discusses  this  procedure  [8]. 


V.  THE  CHARACTERISTIC  SURFACES 


The  characteristic  surfaces  of  the  hyperbolic  diffusion 
equation  provide  a  useful  parameter  of  the  effectiveness 
of  diffusion.  The  characteristics  of  the  hyperbolic 
diffusion  equation  are  the  loci  of  points  of  discontinuity 
of  the  atmospheric  contaminant.  All  solutions  exist 
within  the  characteristics;  there  is  no  concentration  out¬ 
side  of  the  characteristics.  The  characteristics  in 
effect  represent  wave  fronts  of  the  diffusing  particles. 

Pasquill  [25]  finds  that  the  vertical  spread  of  diffus¬ 
ing  particles  in  a  parabolic  Lagrangian  similarity  treatment 
is  best  represented  by  an  extreme  height  encompassing  most 
of  the  particles,  essentially  as  In  the  hyperbolic  treat¬ 
ment  of  Nonin.  In  the  present  hyperbolic  model  this 
concept  is  intrinsic  rather  than  arbitrary  as  in  the 
parabolic  equation  application.  The  extremum  corresponds 
to  the  characteristics  of  the  hyperbolic  equation. 


VI.  THE  PARABOLIC  DIFFUSION  EQUATION  AS  A  LIMITING  FORM 
OF  THE  HYPERBOLIC  EQUATION 

Using  a  method  similar  to  Davies  [8LC93.  and  Monln 
[18],  the  parabolic  diffusion  equation  will  be  derived 
from  the  hyperbolic  diffusion  equation.  The  momentum 
equation  may  be  written 


Hb  ^ 

{ (pV+S).+  $*  CeW+SVeVS+e<vv>]+fa(f ‘Sj+Jd  (J*S)+ky(k*S)-«p  }• 
*  »  0 

o  *33  o  -  0  (23) 

0  0  *kk 

T 


If  1/a,  1/5 ,  1/y»  are  considered  as  parameters,  and  the 
following  limits  are  said  to  exist, 

£il  0  0 

{[i+jeic]*<vv>j*  o  ^33  o  xinmJ-  ♦  ck323  ♦  tK]3k  (24) 

0  0  lick 

lU  0  0 

e 

{♦•«$$>}•  0  *33  0  "v£  +  vj+vk  (25) 

•  ..  11m  1  2  3 

0  0  jkk 


.  & 


ill  0  0 

Cl 

i*  o  ijj  o  "vi4.yj4.vk 

r*  __  llm  «i  *  *2  *3 

0  0  JLkk 

Y 


limits: 


llm 


llm  a, «,Y 
llm  *.<tv 
llm  I-m» 


(26) 
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•  x+&rQ&*y 


then*  taking  the  Units  and  substituting  the  monentun 
equation  into  the  continuity  equation,  a  parabolic 
diffusion  equation  is  found. 


3=  [%<?>»]  ♦ 


J.fckV: 


<f)it  3 


»r. 


where  VDi  -  V±  -  vA  ♦  vel 


(27) 


Essentially  we  have  the  hyperbolic  equation  going  to  the 
parabolic  equation  after  times  very  large  compared  to  - , 

1.  The  parabolic  diffusion  equation  found  above  lsa 
the^three  dimensional  Fokker-Planck  equation  [29]  [30]. 

If  V  -  V(z)i  and  <f$>  has  no  d^agpna].  terms,  equivalent 
to  partition  of  energy  In  the  1,  Jp  k,  directions,  then 
we  have  the  seml-emperlcal  parabolic  diffusion  equation 
for  the  lower  atmosphere  [20]  [24],  If  In  addition  %  *  0, 
the  Smoluchowskl  equation  is  obtained  p]. 


VII 


FUTURE  DEVELOPMENT 


( t  In  the  case  of  constant  mean  velocity,  and  isotropic 
«vv»,  the  present  model  reduces  to  those  of  Davies  [8], 
Goldstein  [13],  and  Monln  [18]  in  the  appropriate  number 
of  spatial  dimensions.  It  is  not  fully  clear  under  what 
conditions  the  models  are  consistent  with  the  various 
spectra.  In  addition  the  effects  of  meandering  have  not 
been  included.  Further  refinement  and  development  of  the 
present  model,  especially  with  respect  to  spectra  and 
pair  particle  density  consideration,  may  be  fruitful  along 
the  phase  space  approach  of  Tchen  [29 ]  and  the  statistical 
representation  of  meandering  by  Gifford  [11]  [12].  These 
approaches  are  consistent  with  the  development  of  the 
present  model.  The  role  of  large  and  small  eddy  scales 
may  then  be  explicitly  Included. 
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